{
    "cells": [
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "### Shared routines\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "import numpy as np\n",
                "import scipy.spatial\n",
                "import matplotlib\n",
                "import matplotlib.pyplot as plt\n",
                "########################################################################\n",
                "#########  Data Generating Functions ###################################\n",
                "########################################################################\n",
                "def generate_sensors(k = 7, d = 2):\n",
                "   \"\"\"\n",
                "   Generate sensor locations.\n",
                "   Input:\n",
                "   k: The number of sensors.\n",
                "   d: The spatial dimension.\n",
                "   Output:\n",
                "   sensor_loc: k * d numpy array.\n",
                "   \"\"\"\n",
                "   sensor_loc = 100*np.random.randn(k,d)\n",
                "   return sensor_loc\n",
                "\n",
                "def generate_data(sensor_loc, k = 7, d = 2,\n",
                "                     n = 1, original_dist = True, sigma_s = 100):\n",
                "   \"\"\"\n",
                "   Generate the locations of n points and distance measurements.\n",
                "\n",
                "   Input:\n",
                "   sensor_loc: k * d numpy array. Location of sensor.\n",
                "   k: The number of sensors.\n",
                "   d: The spatial dimension.\n",
                "   n: The number of points.\n",
                "   original_dist: Whether the data are generated from the original\n",
                "   distribution.\n",
                "   sigma_s: the standard deviation of the distribution\n",
                "   that generate each object location.\n",
                "\n",
                "   Output:\n",
                "   obj_loc: n * d numpy array. The location of the n objects.\n",
                "   distance: n * k numpy array. The distance between object and\n",
                "   the k sensors.\n",
                "   \"\"\"\n",
                "   assert k, d == sensor_loc.shape\n",
                "   obj_loc = sigma_s*np.random.randn(n, d)\n",
                "   if not original_dist:\n",
                "        obj_loc = sigma_s*np.random.randn(n, d)+([300,300])\n",
                "   distance = scipy.spatial.distance.cdist(obj_loc,\n",
                "                                            sensor_loc,\n",
                "                                            metric='euclidean')\n",
                "   distance += np.random.randn(n, k)\n",
                "   return obj_loc, distance\n",
                "\n",
                "def generate_data_given_location(sensor_loc, obj_loc, k = 7, d = 2):\n",
                "   \"\"\"\n",
                "   Generate the distance measurements given location of a single object and sensor.\n",
                "\n",
                "   Input:\n",
                "   obj_loc: 1 * d numpy array. Location of object\n",
                "   sensor_loc: k * d numpy array. Location of sensor.\n",
                "   k: The number of sensors.\n",
                "   d: The spatial dimension.\n",
                "\n",
                "   Output:\n",
                "   distance: 1 * k numpy array. The distance between object and\n",
                "   the k sensors.\n",
                "   \"\"\"\n",
                "   assert k, d == sensor_loc.shape\n",
                "\n",
                "   distance = scipy.spatial.distance.cdist(obj_loc,\n",
                "                                            sensor_loc,\n",
                "                                            metric='euclidean')\n",
                "   distance += np.random.randn(1, k)\n",
                "   return obj_loc, distance\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "## Part B\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "########################################################################\n",
                "#########  Gradient Computing and MLE ###################################\n",
                "########################################################################\n",
                "def compute_gradient_of_likelihood(single_obj_loc, sensor_loc,\n",
                "                                    single_distance):\n",
                "    \"\"\"\n",
                "    Compute the gradient of the loglikelihood function for part a.\n",
                "\n",
                "    Input:\n",
                "    single_obj_loc: 1 * d numpy array.\n",
                "    Location of the single object.\n",
                "\n",
                "    sensor_loc: k * d numpy array.\n",
                "    Location of sensor.\n",
                "\n",
                "    single_distance: k dimensional numpy array.\n",
                "    Observed distance of the object.\n",
                "\n",
                "    Output:\n",
                "    grad: d-dimensional numpy array.\n",
                "\n",
                "    \"\"\"\n",
                "    grad = np.zeros_like(single_obj_loc)\n",
                "    #Your code: implement the gradient of loglikelihood\n",
                "    ### start codeA ###\n",
                "\n",
                "    ### end codeA ###\n",
                "\n",
                "    return grad\n",
                "\n",
                "def find_mle_by_grad_descent_part_b(initial_obj_loc,\n",
                "        sensor_loc, single_distance, lr=0.001, num_iters = 10000):\n",
                "    \"\"\"\n",
                "    Compute the gradient of the loglikelihood function for part a.\n",
                "\n",
                "    Input:\n",
                "    initial_obj_loc: 1 * d numpy array.\n",
                "    Initialized Location of the single object.\n",
                "\n",
                "    sensor_loc: k * d numpy array. Location of sensor.\n",
                "\n",
                "    single_distance: k dimensional numpy array.\n",
                "    Observed distance of the object.\n",
                "\n",
                "    Output:\n",
                "    obj_loc: 1 * d numpy array. The mle for the location of the object.\n",
                "\n",
                "    \"\"\"\n",
                "    obj_loc = initial_obj_loc\n",
                "    # Your code: do gradient descent\n",
                "    ### start codeB ###\n",
                "\n",
                "    ### end codeB ###\n",
                "\n",
                "    return obj_loc\n",
                "    ########################################################################\n",
                "    #########  Part B ########################################################\n",
                "    ########################################################################\n",
                "\n",
                "lr = 1.0\n",
                "# Your code: set some appropriate learning rate here\n",
                "### start codeC ###\n",
                "\n",
                "### end codeC ###\n",
                "\n",
                "np.random.seed(0)\n",
                "sensor_loc = generate_sensors()\n",
                "obj_loc, distance = generate_data(sensor_loc)\n",
                "single_distance = distance[0]\n",
                "print('The real object location is')\n",
                "print(obj_loc)\n",
                "# Initialized as [0,0]\n",
                "initial_obj_loc = np.array([[0.,0.]])\n",
                "estimated_obj_loc = find_mle_by_grad_descent_part_b(initial_obj_loc,\n",
                "        sensor_loc, single_distance, lr=lr, num_iters = 10000)\n",
                "print('The estimated object location with zero initialization is')\n",
                "print(estimated_obj_loc)\n",
                "\n",
                "# Random initialization.\n",
                "initial_obj_loc = np.random.randn(1,2)\n",
                "estimated_obj_loc = find_mle_by_grad_descent_part_b(initial_obj_loc,\n",
                "        sensor_loc, single_distance, lr=lr, num_iters = 10000)\n",
                "print('The estimated object location with random initialization is')\n",
                "print(estimated_obj_loc)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "## Part C\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "########################################################################\n",
                "######### Part c #################################################\n",
                "########################################################################\n",
                "def log_likelihood(obj_loc, sensor_loc, distance):\n",
                "  \"\"\"\n",
                "  This function computes the log likelihood (as expressed in Part a).\n",
                "  Input:\n",
                "  obj_loc: shape [1,2]\n",
                "  sensor_loc: shape [7,2]\n",
                "  distance: shape [7]\n",
                "  Output:\n",
                "  The log likelihood function value.\n",
                "  \"\"\"\n",
                "  # Your code: compute the log likelihood\n",
                "  func_value = 0.0\n",
                "  ### start codeD ###\n",
                "\n",
                "  ### end codeD ###\n",
                "\n",
                "  return func_value\n",
                "\n",
                "\n",
                "########################################################################\n",
                "######### Compute the function value at local minimum for all experiments.###\n",
                "########################################################################\n",
                "num_sensors = 7\n",
                "\n",
                "np.random.seed(100)\n",
                "sensor_loc = generate_sensors(k=num_sensors)\n",
                "\n",
                "# num_data_replicates = 10\n",
                "num_gd_replicates = 100\n",
                "\n",
                "obj_locs = [[[i,i]] for i in np.arange(0,1000,100)]\n",
                "\n",
                "func_values = np.zeros((len(obj_locs),10, num_gd_replicates))\n",
                "# record sensor_loc, obj_loc, 100 found minimas\n",
                "minimas = np.zeros((len(obj_locs), 10, num_gd_replicates, 2))\n",
                "true_object_locs = np.zeros((len(obj_locs), 10, 2))\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "for i, obj_loc in enumerate(obj_locs):\n",
                "    for j in range(10):\n",
                "        obj_loc, distance = generate_data_given_location(sensor_loc, obj_loc,\n",
                "                                                         k = num_sensors, d = 2)\n",
                "        true_object_locs[i, j, :] = np.array(obj_loc)\n",
                "        for gd_replicate in range(num_gd_replicates):\n",
                "            initial_obj_loc = np.random.randn(1,2)* (100 * i+1)\n",
                "            obj_loc = find_mle_by_grad_descent_part_b(initial_obj_loc,\n",
                "                       sensor_loc, distance[0], lr=0.1, num_iters = 1000)\n",
                "            minimas[i, j, gd_replicate, :] = np.array(obj_loc)\n",
                "            func_value = log_likelihood(obj_loc, sensor_loc, distance[0])\n",
                "            func_values[i, j, gd_replicate] = func_value\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "########################################################################\n",
                "######### Calculate the things to be plotted. ###\n",
                "########################################################################\n",
                "local_mins = [[np.unique(func_values[i,j].round(decimals=2)) for j in range(10)] for i in range(10)]\n",
                "num_local_min = [[len(local_mins[i][j]) for j in range(10)] for i in range(10)]\n",
                "proportion_global = [[sum(func_values[i,j].round(decimals=2) == min(local_mins[i][j]))*1.0/100 \\\n",
                "                       for j in range(10)] for i in range(10)]\n",
                "\n",
                "\n",
                "num_local_min = np.array(num_local_min)\n",
                "num_local_min = np.mean(num_local_min, axis = 1)\n",
                "\n",
                "proportion_global = np.array(proportion_global)\n",
                "proportion_global = np.mean(proportion_global, axis = 1)\n",
                "\n",
                "########################################################################\n",
                "######### Plots. #######################################################\n",
                "########################################################################\n",
                "fig, axes = plt.subplots(figsize=(8,6), nrows=2, ncols=1)\n",
                "fig.tight_layout()\n",
                "plt.subplot(211)\n",
                "\n",
                "plt.plot(np.arange(0,1000,100), num_local_min)\n",
                "plt.title('Number of local minimum found by 100 gradient descents.')\n",
                "plt.xlabel('Object Location')\n",
                "plt.ylabel('Number')\n",
                "#plt.savefig('num_obj.png')\n",
                "# Proportion of gradient descents that find the local minimum of minimum value.\n",
                "\n",
                "plt.subplot(212)\n",
                "plt.plot(np.arange(0,1000,100), proportion_global)\n",
                "plt.title('Proportion of GD that finds the global minimum among 100 gradient descents.')\n",
                "plt.xlabel('Object Location')\n",
                "plt.ylabel('Proportion')\n",
                "fig.tight_layout()\n",
                "plt.savefig('prop_obj.png')\n",
                "\n",
                "########################################################################\n",
                "######### Plots of contours. ###########################################\n",
                "########################################################################\n",
                "np.random.seed(0)\n",
                "# sensor_loc = np.random.randn(7,2) * 10\n",
                "x = np.arange(-10.0, 10.0, 0.1)\n",
                "y = np.arange(-10.0, 10.0, 0.1)\n",
                "X, Y = np.meshgrid(x, y)\n",
                "obj_loc = [[0,0]]\n",
                "obj_loc, distance = generate_data_given_location(sensor_loc,\n",
                "                                                 obj_loc, k = num_sensors, d = 2)\n",
                "\n",
                "Z =  np.array([[log_likelihood((X[i,j],Y[i,j]),\n",
                "                               sensor_loc, distance[0]) for j in range(len(X))] \\\n",
                "               for i in range(len(X))])\n",
                "\n",
                "\n",
                "plt.figure(figsize=(10,4))\n",
                "plt.subplot(121)\n",
                "CS = plt.contour(X, Y, Z)\n",
                "plt.clabel(CS, inline=1, fontsize=10)\n",
                "plt.title('With object at (0,0)')\n",
                "#plt.show()\n",
                "\n",
                "np.random.seed(0)\n",
                "# sensor_loc = np.random.randn(7,2) * 10\n",
                "x = np.arange(-400,400, 4)\n",
                "y = np.arange(-400,400, 4)\n",
                "X, Y = np.meshgrid(x, y)\n",
                "obj_loc = [[200,200]]\n",
                "obj_loc, distance = generate_data_given_location(sensor_loc,\n",
                "                                                 obj_loc, k = num_sensors, d = 2)\n",
                "\n",
                "Z =  np.array([[log_likelihood((X[i,j],Y[i,j]),\n",
                "                               sensor_loc, distance[0]) for j in range(len(X))] \\\n",
                "               for i in range(len(X))])\n",
                "\n",
                "\n",
                "# Create a simple contour plot with labels using default colors.  The\n",
                "# inline argument to clabel will control whether the labels are draw\n",
                "# over the line segments of the contour, removing the lines beneath\n",
                "# the label\n",
                "#plt.figure()\n",
                "plt.subplot(122)\n",
                "CS = plt.contour(X, Y, Z)\n",
                "plt.clabel(CS, inline=1, fontsize=10)\n",
                "plt.title('With object at (200,200)')\n",
                "#plt.show()\n",
                "plt.savefig('likelihood_landscape.png')\n",
                "\n",
                "\n",
                "########################################################################\n",
                "######### Plots of Found local minimas. ###########################################\n",
                "########################################################################\n",
                "#sensor_loc\n",
                "#minimas = np.zeros((len(obj_locs), 10, num_gd_replicates, 2))\n",
                "#true_object_locs = np.zeros((len(obj_locs), 10, 2))\n",
                "object_loc_i = 5\n",
                "trail = 0\n",
                "\n",
                "plt.figure()\n",
                "plt.plot(sensor_loc[:, 0], sensor_loc[:, 1], 'r+', label=\"sensors\")\n",
                "plt.plot(minimas[object_loc_i, trail, :, 0], minimas[object_loc_i, trail, :, 1], 'g.', label=\"minimas\")\n",
                "plt.plot(true_object_locs[object_loc_i, trail, 0], true_object_locs[object_loc_i, trail, 1], 'b*', label=\"object\")\n",
                "plt.title('object at location (%d, %d), gradient descent recovered locations' % (object_loc_i*100, object_loc_i*100))\n",
                "plt.legend()\n",
                "plt.savefig('2D_vis.png')\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "1+3\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": []
        }
    ],
    "metadata": {
        "kernelspec": {
            "display_name": "Python 3",
            "language": "python",
            "name": "python3"
        },
        "language_info": {
            "codemirror_mode": {
                "name": "ipython",
                "version": 3
            },
            "file_extension": ".py",
            "mimetype": "text/x-python",
            "name": "python",
            "nbconvert_exporter": "python",
            "pygments_lexer": "ipython3",
            "version": "3.8.3"
        }
    },
    "nbformat": 4,
    "nbformat_minor": 4
}